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ABSTRACT: We present two continuum models A and B to study the con- 
vective instability of granular materials subjected to vibrations. We carry out 
the linear stability analysis for model A and uncover the instability mecha- 
nism as a supercritical bifurcation of a bouncing solution. We also explicitly 
determine the onset of convection as a function of control parameters. The 
simulations results are in excellent agreement with the stability analysis. 
Additional feature of the model B is the inclusion of the relaxation term in 
the momentum equation, which appears to be crucial in capturing what is 
missing in model A, in particular, in reproducing experimental convection 
patterns for large aspect ratio, both horizontally, in which case convective 
rolls move toward the surface, and vertically in which case convective rolls 
survive near the wall but are suppressed in the bulk region. 



1. INTRODUCTION 

It was Faraday who discovered the convective instability in a vibrated 
granular bed in 1831. Initially the flat surface of the granular pile develops 
a heap upon vibrations, along the surface of which grains roll down causing 
small or large scale avalanches. Once formed, such a heap is stable, because 
of the simultaneous formation of permanent convective rolls inside the heap. 
Unlike Rayleigh-Bernard convection in fluids, however, the origin of this in- 
stability has remained relatively unexplored since its discovery, but recently 
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two simultaneous push from experimental side ( Clement et al, 1992; Pak and 
Behringer, 1994) with the use of MRI or X-ray method (Knight et al, 1993) 
from large scale computer simulations based on the distinct element method 
(Taguchi, 1992; Gallas et al, 1992) have aided our understanding through 
visualization. However, the theoretical efforts(Haff,1983; Bourzutschky and 
Miller, 1995) to uncover the basic mechanism of this convective instabil- 
ity have not been remarkable, still largely focused on producing convective 
patterns through computer models and simulations. We have recently under- 
taken steps, based on two continuum models, to remedy this situation, which 
appear to have captured the essence of granular convection. Considering a 
potentially important industrial application of size segregation and a recent 
evidence (Knight et al, 1993) of the convection connection in conjunction 
with the conventionally held reorganization of grains, we consider the search 
for the origin of granular convection quite important. This is a brief sum- 
mary of our effort along this direction. For details, see Hayakawa et al(1995) 
and Yue(1995). 



2. MODEL A 

We have studied two models. Both are based on Navier-Stokes type 
continuum eqs, but ignore temperature equation assuming the existence of a 
global temperature throughout the bed. We first present the model A. 

The staring point of model A is the recognition that the most fundamental 
aspect of the vibrating bed, apart from the obvious fixed bed solution with 
no external driving, is the existence of a uniform bouncing of a collection of 
particles, a solid or a fluidized block with no internal degrees of freedom. In 
such a case, the bouncing solution, represented by the motion of a ball on 
a vibrating platform, satisfies z = ( — 1 + T sin t)9( — 1 + Tsini) where z is 
the vertical coordinate of the granular block and the 8(x) = 1 for x > and 
6(x) = for otherwise. Next, in the presence of internal degrees of freedom 
such as rotation and/or translation, we define two coarse-grained dynamical 
variables: the density p(r, t) and the velocity v(r,t) of the granular system. 
In the box fixed frame, eq. then modifies into: 

d t p + V-(pv) = (1) 
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<9 t v + (v-V)v = z(rsin t — l — X) VP 

P 

+^[V 2 v + X V(V-v)] (2) 

where z is the unit vector in the vertical direction and A is the Lagrange 
multiplier. A = for free motion and A = Tsint — 1 for stationary state. 
Note that the first term in the right hand side of (@) is due to the uniform 
bouncing and the third term is the energy dissipation effectively represented 
by the Reynolds number R and the bulk viscosity x- The pressure term P 
requires some discussion (Hayakawa et al, 1995) but the Van der Waals model 
P = Tp/(1 — bp) is a reasonable choice, where T represents the effective 
temperature which might be a global variable and b is a constant of order 
unity. 

To check the validity of our picture, we have solved ([!]) and (g) numeri- 
cally in two dimension with no slip boundary conditions at the side walls as 
well as at the top and the bottom plates. Note that the top plate suppresses 
complicated surface motion of vibrating beds and allows us to use the simpli- 
fied picture. Since the granular fluid is confined in a box, we do not introduce 
A explicitly in the simulations. The absence of A and the presence of the top 
wall is expected to cause the appearance of the bouncing solution for T < 1 
in contrast with the real situation but its omission would not change the 
essence of the dynamics. In the same spirit, we have ignored x an d b in our 
simulations. 

For T < T c , the bouncing solution is expected to appear inside the bed 
and the density and the velocity at a given point oscillates with the same 
frequency of the vibration. ( Fig. 2) Upon increasing T further to T = 1.2, 
which is beyond the predicted T c = 1.12 determined by (Jm(Xc) = 0, we find 
that the bouncing solution has disappeared and the permanent convective 
rolls have developed inside the bulk (Fig.3). The wavelength of the most 
unstable mode by the linear stability analysis is about q m ~ 0.4, which is not 
far from the actual wavelength of the convective rolls: q = 2n / 'X = 2ir / 'L w 
0.6. 



3. MODEL B 
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We now introduce another model as Model B. Although the mass conser- 
vation remains same as in ([[]), the momentum conservation (Q) now changes 
into 

d t v x + {v-V)v x = -{cl/ p)d xP + iiV 2 v x (3) 
d t v z + (v-V)v z = (V(p,t)-v z ))/r 

-{cl/p)d zP + pV 2 v z (4) 

where Cq ~ T is the sound speed. The difference between model A and B is 
the presence of a relaxation term in the z direction (|4]), which is represented 
by an average function V(p) with the relaxation time r. The origin of such 
a term has been discussed in Hong et al(1994) in an attempt to introduce 
correlations among grains or voids in the diffusing void model(DVM). In 
the DVM, the void speed is only a function of the local density, namely 
v z = V(p) + diffusion term. However, a void is a compressible hydrodynamic 
object that changes and adjusts its shape to conform to the surrounding, not 
instantaneously, but in a given time. So, it may be more appropriate to write 
down the time dependent equation for the velocity in a manner given by (^) 
than simply assuming a fixed value at a given local density. The presence 
of such relaxation process may be effectively equivalent to assuming a drag 
force acting on a void. 

These coupled equations (H) and (|) are also known as the traffic model or 
two-phase model for fluidized beds that have been widely used for mixtures 
of gas and granular particles. Functions in the model may be inferred from 
the Enskoq equation; namely —v z /r is the drag term imposed externally on 
the particle. In the case of no interstitial fluid, its origin lies in the frictions 
of the front and the rear glass of the container and from the wall. Further, 
the Enskoq pressure, Tp(l + f(p)p/2) with f(p) the correlation function, 
produces an extra term V(p) in addition to the hard sphere pressure T p. In 
this case, the coefficient of V(p) is proportional to the gravity g. The net 
effect is for the void (or particle) to adjust its speed, v z , around the average 
value V(p) in a given time r. While deriving the exact form for the function 
V(p) is nontrivial, we know it must be a decreasing function of density and 
have a cut off at the closed packed density p c . Hence, we have chosen a 
simple form: V(p) = V r D (p)(— 1 + rsin(u;t)) V (p) = (p c — p)^9(p c — p) with 
9 function and /3 < 1. We have investigated the eqs.(^) and (^) numerically. 
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The mechanism for the convection seems to be similar to model A, namely 
the superciritcal bifurcation of a bouncing solution. However, two notable 
differences emerge. First, when the aspect ratio increases vertically from, 
say one to two, the convection rolls that initially occupied the whole box 
move toward the surface as shown in Fig. 3 and the motion of the particles 
are fairly confined near the surface. This is consistent with the experiments 
and MD simulations. 

Second, when the aspect ratio increases horizontally, then convective rolls 
inside the bulk are suppressed and they appear only near the wall, which is 
not shown here. This again is consistent with the MD simulation results of 
Taguch(1992). Hence, the role of drag appears to be important in granular 
convection. 

4. DISCUSSION 

First, the bouncing solution as a basic state for granular convection seems 
to have been confirmed in the simulations of both models A and B. 

Second, the role of boundary conditions. We have employed no slip 
boundary conditions at the walls and the plate. Further, we have put the 
rigid wall at the top and thus suppressed the surface motion. Granular mate- 
rials have been shown to exhibit very different motion near the wall in a zet 
flow experiment under gravity(Caram and Hong, 1993) and there has been 
some attempt to use the negative or positive slip to control the convective 
patterns (Bourzutschky and Miller, 1995). More detailed studies to derive 
reasonable boundary conditions at the wall are required. 

Third, the role of interstitial fluid. Our model A assumes no interstitial 
fluid such as air, and it predicts a series of rolls for the vertically large aspect 
ratio. Model B on the other hand predicts that the convection is suppressed 
in the bulk region but is confined near the surface, which is in accordance with 
experiments and with the results of the two-fluid model with an interstitial 
fluid. Perhaps, the origin of drag, whether it is coming from friction at the 
walls, or from the viscous effect of the interstitial fluid, may not be relevant 
once it is present. The suppression of the convection in the bulk is due to 
the locking mechanism of grains for near closed packed density, which was 
taken into account in model B by a cut off in V(p), namely V(p) =0 for 
P > Pc- We need more extensive studies of both models A and B to make 
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quantitative comparison with experiments. 

HH is, in part, supported by Grant-in- Aid for Scientific Research from 
Japanese Ministry of Education, Science and Culture (No. 08226206 and 
08740308). 

REFERENCES 

Bourzutschky, M. & Miller, J. 1995. Granular convection in a vibrating fluid. 
Phys. Rev. Lett. 

74: 2216-2219. 

Caram,H., & Hong, D. 1992. Diffusing void model for granular flow. Mod. 
Phys. Lett. B. Vol.(6):761-771. 

Clement, E., Duran J., & Rajchenbach,J. 1992. Experimental study of heaping 
in a two dimen- sional sand pile. Phys. Rev. Lett. 69:1189- 1192. 

Gallas,J., Hermann, H.,& Sokolowski,S. 1992. Convection cells in vibrating 
granular media. Phys. Rev. Lett. 69: 1371-1374. 

Haff,P. 1983. Granular flow as a fluid mechanical phenomena. J. Fluid. 
Mech. 134: 401 

Hayakawa,H,S. Yue,& Hong,D. 1995. Hydrodynamic description of granular 
convection. Phys. Rev. Lett. 75:2328-2331. 

Hong et al. 1994. Granular relaxation under tapping and the traffic equation. 
Phys. Rev. E. 50: 4123-4135. 

Knight, J., Jaeger. H, & Nagel,S. 1993. Vibration induced size separation in 
granular media: convection connection Phys. Rev. Lett. 70: 3728-3731. 

Taguch,Y-h. 1992. New origin of a convective motion. Phys. Rev. Lett. 69: 
1367-1369. 

Yue,S 1995. Nonlinear phenomena in materials failure and granular dynam- 
ics. Ph.D thesis, Lehigh University. 



6 



Figure Captions 



Figure l:The effective growth rate a e ff(q) as a function of the wave number 
q for T = 1.05(diamond),for which a e ff(q) < for all values of q, while for 
T = 1.2 > T c = 1.12, cr e ff(q) becomes positive for a band of q(square). r c 
is determined by the condition that the maximum of cr e ff(q) becomes zero 
at T c . (cross) The parameters used are: T e = R = 10 and L = 10. (Figure is 
missing. Please see PRL 75,2328, 1995) 

Figure 2: For T = 1.2 > T c = 1.12, the bouncing solution becomes unstable 
and the permanent convective rolls appear inside the box. The arrows are 
the velocity vectors pointing upward. 

Figure 3:For T = 1.2 > T c and for large aspect ratio along the vertical, the 
convective rolls move toward the surface. 

Figure 4:For T = 1.2 > r c and for large aspect ratio along the horizontal, the 
convective rolls in the bulk are suppressed and survive only near the wall. 
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